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A Bose-Hubbard model, describing bosons in a harmonic trap with a superimposed optical lattice, 
is studied using a fast and accurate variational technique (MF+NRG): the Gutzwiller mean-field 
(MF) ansatz is combined with a Numerical Renormalization Group (NRG) procedure in order 
to improve on both. Results are presented for one, two and three dimensions, with particular 
attention to the experimentally accessible momentum distribution and possible satellite peaks in this 
distribution. In one dimension, a comparison is made with exact results obtained using Stochastich 
Series Expansion. 
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I. INTRODUCTION 

The recent experiments by Greiner et al. pj on bosons 
confined in an optical lattice demonstrated the transition 
between a Mott phase and a Superfluid phase (SF), as 
was first predicted by Jaksch et al. 0j- The experiments 
are adequately described by a single band Bose-Hubbard 
Hamiltonian with on-site repulsion only. This type of 
repulsion leads to two different phases. A Mott insu- 
lating phase can exist at commensurate fillings, with a 
quantum phase transition to a superfluid as the density 
is shifted or the interaction strength weakened. In the 
experiments however, the quadratic confining potential 
adds a new term to the Hamiltonian that cuts off any 
long range correlations, but Mott and superfluid regions 
can still occur. The experiments led to a complete revival 
of interest in the bosonic model, thanks to the unprece- 
dented control over the physical parameters compared to 
former realizations. 

Other experimental realizations of bosonic lattice sys- 
tems include 4 He on graphite [j| , superconducting islands 
or grains connected by Josephson junctions Q. In this 
case, Cooper paired fermions are considered as bosons, 
at least approximately. Recently, attempts have been 
made to investigate at what scales the fermionic nature of 
paired fermions can play a role p| , and the result is that 
for the energy scales considered here, individual atoms 
can safely be described by a one-boson operator b>, as 
expected. 

Since the work by Fisher et al. [||, the determina- 
tion of the ground state phase diagram of the Bose- 
Hubbard model with on-site repulsion only has at- 
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tracted a lot of attention. Analytic studies using mean- 
field theory 0, and renormalization group tech- 
niques @ led to a deeper physical understanding of the 
model. Strong coupling expansions 0, gave a better 
quantitative pict ure, while quantum Monte-Carlo sim- 
ulations 0, Il2l Il3l fl4L Il5| were carried out in one 
and two dimensions. The one-dimensional case was re- 
cently investigated using the Density Matrix Renormal- 
ization Group (DMRG) 0, 0, yielding the most ac- 
curate results at present. Longer range interactions can 
cause c harg e density wave, stripe or even supersolid or- 
der 0, Il8l . Disorder and impurities can have dra- 
matic effects la, 

m mm Hum and lead to even other 

phases. The model with a quadratic confining poten- 
tial 2] has been addressed in one dimension [23] and for 
a small lattice in three dimensions |24j . using quantum 
Monte-Carlo methods. 

In view of the enormous success of DMRG |25| in 
bosonic 01 an d fermionic [2(j real-space lattice mod- 
els in one dimension, DMRG has been extended beyond 
these models, towards app lications in metallic grains |27j . 
quantum chemistry |28t I29I ] and first attempts have 
even been undertaken towards applications in nuclear 
physics [3(J. Unfortunately, an exact DMRG study in 
all 3 dimensions of the Bose-Hubbard model is not fea- 
sible with current computer power. DMRG was a sub- 
stantial improvement |25| on the older Numerical Renor- 
malization Group (NRG), which had only a poor rep- 
utation in dealing with long-range interactions between 
fermions |3l|, 1321 . but we found it useful for bosonic sys- 
tems (see also |33|). 

The basic philosophy of this work consists of extend- 
ing one-site mean-field theory to larger blocks, following 
the idea of the Renormalization Group method. A fully 
variational method is obtained which incorporates corre- 
lations beyond mean-field at low computational cost and 
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which is able, unlike NRG, to accurately describe the 
different phases of the Mott-SF transition. We judged 
the computational cost as primordial, so that extensions 
to large lattices are in reach and so that a direct simula- 
tion of the experimental parameters can be accomplished, 
while the computational uncertainties remain well within 
the experimental uncertainty range. A disadvantage of 
the method is the breaking of number conservation dur- 
ing the intermediate steps of the renormalization proce- 
dure. In the final step particle number should be restored 
in principle, but this restoration is only partial when an 
insufficient amount of states are kept during the model 
space truncation. The quantity that directly relates to 
experiment is the momentum distribution. It was argued 
in Ref. .24] that the appearance of satellite peaks in the 
momentum distribution signals the appearance of a Mott 
region in the center of the trap. We will especially focus 
on the momentum distribution and on the issue of the 
satellite peaks. 

The organization of the paper is as follows: In Sec. ITT1 
we explain the method, in Sec. II I II compare it to exact 
diagonalization methods for small latices and in Sec. IIVI 
we present results in 1, 2 and 3 dimensions. We end with 
the conclusion and the acknowledgments. 



II. A NUMERICAL RENORMALIZATION 
GROUP METHOD 

Bosons in an optical lattice realize a Bose-Hubbard 
Hamiltonian and more specifically, we consider here 
the soft-core Bose-Hubbard Hamiltonian in the grand- 
canonical ensemble in d dimensions and subject to a con- 
fining field, 

H = -tJ2 b\b, + \lJ £ rn(ni - 1) - jiK (1) 

(id) i » 

The sum is over the nearest neighbors only. The 

operator b\ creates a boson on lattice site i while hi re- 
moves it. The operator counts the local density on 
site i. The operator N denotes the total number opera- 
tor, N = n i- We take the distance a between adjacent 
sites equal to a = 1. The confining field acts as a local 
site dependent chemical potential that can be added 
to the chemical potential /i to form /ij = ej — fi. In case 
of a site dependent \Xi we speak of the inhomogeneous or 
confined model, in case of a uniform /i we speak of the 
unconfined or homogeneous model. We also define the 
coordination number z = 2d as the number of neighbors 
of each site and consider a linear, square or cubic lattice 
of length L along each axis. The energy scale is set by 
setting t — 1. 

This Hamiltonian is the easiest bosonic model in which 
two different effects compete: the kinetic energy is diago- 
nal in momentum space and tries to delocalize the parti- 
cles over the sites, while the potential energy is diagonal 
in coordinate space and localizes the particles. 
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FIG. 1: The mean-field phase diagram in the units we adopt 
throughout the paper. The Mott lobes are indicated and sur- 
rounded by a superfluid. 

We first discuss the physics of the model in absence 
of disorder in one dimension [|| 0| • When the potential 
energy dominates, the system forms a Mott insulating 
phase at integer densities, that remain pinned at these 
integer values and the phase is incompressible. The Mott 
lobes are surrounded by compressible SF phases, where 
the densities fluctuate. This is visualised in the mean- 
field phase diagram Fig. which can easily be calcu- 
lated @, 13 ■ Note that this phase diagram is approx- 
imate, in the true phase diagram the n = 1 lobe e.g. 
extends to smaller values of U and the lobe closes in a 
point-like fashion. There are two different phase transi- 
tions possible. When keeping the density constant at an 
integer value, phase fluctuations dominate and the tran- 
sition is of the Berezinskii-Kosterlitz-Thouless (BKT) 
type. This transition can only occur at the tip of the 
insulator lobe, that as a consequence closes in the point- 
like fashion. The generic phase transition is driven by 
density fluctuations and belongs to a different universal- 
ity class. For a general d dimension, the BKT transition 
generalizes to the (rf + 1) dimensional XY universality 
class. 

A first approximation for the Bose-Hubbard model is 
the Gutzwiller variational ansatz 0, leading to a de- 
coupling of the individual lattice sites and to a mean-field 
theory (MF). This assumption can be described by the 
following substitution: 

bibj^wi+tfbj-w:, (2) 

with = (pi). This leads to a model where particle 
number symmetry can be broken, and that can exhibit a 
superfluid and a Mott-insulator behavior. 

In order to obtain a higher accuracy, we extend the 
MF approximation to a a Renormalization Group proce- 
dure by taking more correlations into account. It works 
as follows. Just as in MF, first break down the entire 
lattice to single sites and solve the problem for each site 
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separately. The Hilbert space is truncated so that only 
a few basis states are kept on every site. In NRG this 
state selection is based on energy solely, meaning that we 
keep the N s eigenstates corresponding to the N s lowest 
energy eigenvalues. The two sites are combined now to 
form a small block. At this stage, the MF approximation 
( |2J for the hopping term between the two sites can be 
canceled by adding a term (tp* — b\)(ipj —bj), after which 
the Hamiltonian for the two sites becomes 



H 12 =H 1+ H 2 + V (V>*-foi)(^-&,)+h.c (3) 



E 



Here, Hi and H 2 denote the Hamiltonians of the left and 
right site, respectively, and the sum runs over adjacent 
sites that each belong to a different block. In this first 
step, just the two sites 1 and 2 are meant. The Hamil- 
tonian H12 is diagonalized in the space spanned by the 
product states, that are constructed from the individual 
basis states of each site. After diagonalization, only a 
few states are kept again. Physical observables require 
now a rotation, since we have performed a basis rota- 
tion. The procedure repeats itself: the small blocks can 
be joined to form larger blocks which will themselves be 
the building blocks of still larger blocks etc. This proce- 
dure is very similar to DMRG [25| . The main differences 
are that in DMRG the selection of the states is based 
on the eigenvalues of the density matrix instead of on 
the lowest energy values, and secondly that in NRG one 
combines blocks (exploiting symmetry), while in DMRG 
one extends the blocks site by site. In NRG one performs 
one calculation till the lattice is entirely built up, while 
in DMRG one sweeps again through the lattice till con- 
vergence is obtained. DMRG yields results with a higher 
accuracy, but its computational time and memory cost 
requirements are beyond current computer power for di- 
mensions higher than one. 

A new idea is that we improve on the standard NRG 
procedure by adding source terms to the Hamiltonian on 
the edges of the blocks. These terms compensate for the 
interaction with the other blocks in a mean-field way. 
In this way, the Hamiltonian of the local block feels al- 
ready an average contribution of the blocks that have 
not yet accounted for, and that have a non-local influ- 
ence on the block under consideration. After the two 
joining blocks are taken together, these terms need to be 
extracted again. If it were possible to work in an infinite 
Hilbert space the net effect of these source terms would be 
zero, but in a truncated space the calculation will depend 
on the values of these terms. E.g. suppose we are looking 
for a Mott phase and these source terms are set to finite 
values, then we will not find the Mott phase if the source 
term yields contributions to states higher than the cut- 
off. This surely will be the case near the boundary of the 
Mott lobe in a homogeneous system. We have also tried 
to apply improved periodic boundary conditions by use of 
such source terms, contrary to DMRG where one usually 
adopts open boundary conditions. (In general, periodic 
boundary conditions are easier for finite-size scaling.) We 



will come back to the issue of source terms in the next 
section. 

Some operators, such as the total number operator 
squared N 2 also acquire a contribution from the cross- 
terms between the two building blocks. So, more than a 
simple rotation is needed in this case, and contributions 
from the total number operator TV must be taken into ac- 
count. Schematically, (Nf 2 ) = (N?) + (JV|) + 2(N 1 N 2 ), 
in which the indices 1, 2 indicate the two joining blocks. 

One can obtain the one-body density matrix in coor- 
dinate space with the MF+NRG method, but at a low 
accuracy because correlations b\bj are inaccurate when i 
is not the first site of the renormalization procedure, and 
we need the entire matrix for a confined system. How- 
ever, we are able to directly calculate the diagonal of the 
momentum density operator pfc,fc = Pk (in one dimen- 
sion), 



Pk 



hJ 



= co8(k(i-j))btb:j. 



(4) 



Here we sum over all sites of the blocks, easing the prob- 
lem encountered with the one-body density matrix in co- 
ordinate space. This operator can be rewritten as 

= \_. cos(fcz) co$(kj)b\bj + sin(fci) sm(kj)b\bj 



^2cos(ki)b\ ^cos(fcj)6j 



\ieL 



+ (5>n(fcz)6M (^ S in(fei)6 j 

\ieL I \j£R 

Cfcr C*H + Skr Skit- 



(5) 



Hence we have to keep track of linear combinations of 
the creation c] and annihilation operators. The parts 
in which both sites i and j belong to the same block 
(left(L) or right(R)) have been omitted, since their up- 
dating consists only of a rotation to the newly truncated 
basis. When the sites belong to different blocks, there 
is also a contribution of the cross-term just as with the 
operator N 2 , but it still suffices to update the C and S 
operators. The extension to higher dimensions of Eq.© 
is straightforward. We normalize the Fourier transform 
by adding prefactors i- so that the trace of the density 
matrix in momentum space yields the number of particles 
in the system. 



III. COMPARING THE METHOD WITH 
EXACT RESULTS FOR SMALL LATTICES 

In this section we consider a small lattice in one di- 
mension and in the absence of any kind of disorder. 
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TABLE I: Comparison of the energies (-B16 and -B32) per site 
obtained by the MF+NRG method for a modest number of 
states (16 and 32) kept after each diagonalization with the 
results El of a Lanczos diagonalization procedure for a small 
lattice of 8 sites in 1 dimension. The deviations (D16 and 
D32) are indicated and can be made smaller by keeping more 
states after each diagonalization. The mean-field values are 
in the last column. 



u 




E L 


Em 


D 16 (%) 


-E32 


D 32 (%) 


Emf 


2.0 


-0.5 


-1.359 


-1.337 


1.57 


-1.347 


0.91 


-1.24 


4.0 


0.7 


-0.932 


-0.901 


3.38 


-0.914 


1.99 


-0.78 


6.0 


1.8 


-0.656 


-0.612 


6.65 


-0.635 


3.15 


-0.44 


8.0 


2.5 


-0.494 


-0.467 


5.60 


-0.483 


2.40 


-0.27 


10.0 


3.5 


-0.397 


-0.382 


3.53 


-0.392 


1.09 


-0.10 


12.0 


4.0 


-0.331 


-0.324 


2.35 


-0.330 


0.31 


0.00 


14.0 


5.0 


-0.284 


-0.283 


0.50 


-0.284 


0.23 


0.00 


16.0 


6.0 


-0.249 


-0.247 


0.38 


-0.248 


0.07 


0.00 


18.0 


7.0 


-0.221 


-0.221 


0.30 


-0.222 


0.04 


0.00 


20.0 


15.2 


-0.199 


-0.181 


8.70 


-0.199 


0.17 


0.00 



We have checked the code by comparing the resulting 
energies to direct Lanczos diagonalization values for a lat- 
tice containing 8 sites in one dimension. The one dimen- 
sional Bose-Hubbard model with periodic boundaries is 
a worst-case scenario for our MF+NRG procedure. The 
results are summarized in Table [I] The parameters in 
the table vary from a SF phase to a Mott phase. As ex- 
pected, very deep in a Mott phase or in a SF phase we 
obtain a very good accuracy. Note that the Lanczos di- 
agonalization was performed with a fixed boson number, 
while in the MF+NRG we adjusted the chemical poten- 
tial in order to fix the density. The deviations should be 
interpreted accordingly. 

We have also checked observables like the local density 
rii and local compressibility Ki against results obtained 
with the Stochastic Series Quantum Monte-Carlo [35[ 
(SSE) method for larger lattices. Because the calcula- 
tion of the momentum distribution seems most critical, 
we have explicitly shown in Fig.|3the good agreement be- 
tween the calculation of the momentum distribution with 
the renormalization group and the SSE method for a SF 
and a Mott phase. The one-body density matrix with the 
SSE method has been obtained by applying the idea of 
Ref. j3|| to soft-core bosons. The SSE method served as 
a testing ground for the renormalization method here. So 
we have shown that a lot of physics might be examined 
with our new method. 

The parameter N s that fixes how many states are kept 
in the truncation determines the accuracy of the results. 
As we have seen in Table [I] energies can decrease by 
increasing N s , while for observables like the local den- 
sity the fluctuations become smaller. We have examined 
how the grand-canonical potential <& decreases when N s 
is increased in the upper part of Fig. [3] for a system of 
L = 1024 sites (so that finite-size effects can be filtered 
out of this discussion) in the SF phase but very close 
to the generic Mott phase transition. This corresponds 
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FIG. 2: Checking the momentum distribution obtained with 
MF+NRG to a SSE calculation for an unconfined system in 
one dimension of 32 sites. Calculations have been done for 
a system in the Mott phase (U — 6, /1 = 2, dashed line) and 
for a system in the SF phase (U = 2, fi = 1, full line). The 
errors on the SSE data points are shown but very small. In 
the inset, the MF+NRG data points are indicated explicitly 
by small circles, while "+" signs with error bars indicate the 
SSE data points. 



to the worst case scenario for our method. Inclusion of 
just a few states leads to a rapid decrease in the grand- 
canonical potential, but once more than 20 states are 
kept, the potential decreases only very gradually. The ex- 
act result <fr/L = —1.917(2) in Fig.[3]was again obtained 
by the SSE-method, while with N s — 40 we reached 
$/L = —1.90. Without source terms, we found that 
the calculated average grand-canonical potential per site 
was $/L = —1.87 with N s = 40, giving further evidence 
of the usefulness of the source terms. It is the sweeping 
property of the DMRG algorithm that could improve the 
results here substantially, something we tried to avoid 
from the onset since this property is computationally too 
costly in higher dimensions. The discrepancy with the ex- 
act result in Fig. reduces rather slowly at higher values 
of N s , primarily because of the effects of block extension 
(reflected in the curves of the local densities and local 
energies in Fig. 0] for the same effect) and complications 
due to the periodic boundary conditions. The exact re- 
sult should be recovered in the limit of N s equal to the 
dimension of the Fock space for each block. In addition, 
for a system that is already deep in one of both phases, 
the MF+NRG method converges very rapidly to the ex- 
act result, as the energy curve shows in the lower part of 
Fig. G] for a system in the Mott phase. Here, the energy 
and grand potential differ only by a constant. 

On the other hand, the parameter A^ s also largely de- 
termines the required computer time: observables scale 
as Ng per lattice site in memory cost, and the most time 
consuming operation is the rotation of variables, which 
scales as Nf (multiplication of matrices of order N^). All 
our calculations have been performed on a Pentium IV, 
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FIG. 3: Upper: The system is variational in the average 
grand-canonical potential per site 4>/L = ((E) — fj.{N))/ L. As 
the number of kept states N s is increased, the grand-canonical 
potential decreases. The full line is a guide to the eye. The 
parameters are chosen such that the system is near a generic 
phase transition on the SF side (U = 4, fj, = 1, L = 1024). 
This corresponds to the worst case scenario for our method. 
The mean-field (N a ~ 1) and exact (corresponding to N a — 
oo) results are indicated. Lower: The exponential conver- 
gence of the energy is shown. The system is deep in the Mott 
phase with parameters U = 30, [i = 16 and L = 1024. 



1.6GHz or a Pentium III, each with 500MB RAM. Larger 
lattices and higher values of N s can straightforwardly be 
implemented on more performant hardware, but requir- 
ing that all occupation numbers of the truncated states 
are arbitrarily small on the one hand and on the other 
hand wishing to study large lattices in high dimensions 
near a quantum phase transition is still not achievable. 

One of the crucial parameters of the method is the 
source term that is inserted at the boundary of each 
block. If we set it to zero, our method reduces to the 
standard NRG method. It can be seen in Fig.0|that the 
fluctuations can be damped much better in a SF phase 
if we set the source term equal to the MF expectation 
value of the operator a, while the total energy deviates 
now 0.9% from to the exact result instead of 1.4% with- 
out the source terms. For a homogeneous model in the 
thermodynamic limit, the value of (cj) should be site- 
independent and the source terms should be chosen equal 
as well. 

As argued in the previous section, the source terms 
are optional and need to be chosen carefully. It is well 
known [3_^| that a Mott phase can only be found if 
[H,N] = 0. Source terms might violate this condition 
near a Mott-SF transition. The addition of source terms 
might lead to an incorrect prediction of a SF phase when 
the compensation of the source terms in the renormal- 
ization scheme is not complete. The source terms could 
yield contributions to states that are thrown away after 
truncation of the Hilbert space and these contributions 
can be quite large when the parameter N s is chosen too 



FIG. 4: The figure shows how source terms can improve the 
calculations. Local energy per site El, local compressibility 
Ki and local density rn from bottom to top are plotted as a 
function of the lattice index i for a homogeneous model of 8 
sites in a SF phase (U = 2, fi = —0.5, the same values as the 
upper row in Table 0. The dashed line has source terms set 
to zero, while for the full line they are set to their mean-field 
values. 



low. This might lead to an incorrect value of the tran- 
sition point. In addition, even if the transition point of 
a generic Mott-SF phase transition was known exactly, 
and we would study the SF side of this transition, the 
source terms should still be chosen carefully. This can be 
understood as follows: any net contribution of a source 
term will deal with long-range correlations in the same 
way as mean-field does, and we know that the correlators 
predicted by mean-field are only valid in d = 3 dimen- 
sions. E.g. the density is a valid order parameter for the 
generic Mott-SF transition j^l, with 

n ~ /i 1/2 d = 1 

n ~ /i d — 3. (6) 

Briefly said, the source terms would in one and two 
dimensions lead to an improved mean-field theory, in the 
sense that the correlators would approximately have the 
same exponents as in mean-field theory and the Mott lobe 
would extend a little bit farther into parameter space. 
This is also explained in Fig. [3] When N s is large enough, 
these possible dangers become less severe. In the confined 
case, any long-range correlations are effectively cut off 
and the addition of source terms is always expected to 
improve the calculations, as has been verified. 

It was also tempting to study what happens if the 
blocks were extended by a single site only. That leads to 
the same number of diagonalizations but more rotations 
are needed. For the unconfined case, this yielded quite 
good results, often smoother than in the Block Renor- 
malization case. However, for the confined case this pro- 
cedure did not produce regions with integer density and 
should hence only be used with great care. 
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FIG. 5: Evolution of the density in the neighborhood of the 
generic phase transition between the Mott and SF phase for 
one ("+" marks) and three (empty circles) dimensions. The 
dashed line is a fit according to Eq.(|SJl. In these calculations 
the source terms are set to zero and lattices of size L = 1024 
sites were studied. Inclusion of the source terms in one di- 
mension would lead to a similar plot as in the 3D case. 
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FIG. 6: Local density (upper curves) and local compressibil- 
ity (lower curves) for parameters U/z = 5.0 and fjb/z = 2.0 
in one (full line), two (dotted line) and three (dashed line) 
dimensions. MF theory predicts a SF phase, while the true 
phase in one, two and three dimensions is a Mott phase. 



The MF+NRG procedure also offers a substantial im- 
provement over MF results. The MF transition between 
the SF phase and the Mott phase is independent of the 
dimension of the system and is located at U c f» 5.83, 
while a DMRG [rjj study locates it at U c /z rs 3.36 in 
one dimension and a strong-coupling expansion ^(j lo- 
cates it at U c /z « 4.18 in two dimensions. We have 
performed a simulation at U/z — 5.0 in one, two and 
three dimensions. While MF theory predicts a SF phase 
(calculation yields (c) — 0.496) the true phase should be 
a Mott phase in one and two dimensions and we even 
found a Mott phase in three dimensions. These results 
can be seen in Fig. [fj] where the local density and local 
compressibility are shown. 



IV. RESULTS 

Now that we have critically examined the approxima- 
tions made in the MF+NRG, we apply it to parameter 
regions where the results are unambiguous. In all calcula- 
tions the size of the lattices corresponds to the maximal 
achievable size while a sufficient amount of states has 
been kept. 



A. Results in one dimension 

We have, in complete analogy with Ref [23j, confined 
a Bose gas in a lattice of 128 sites in a trapping potential 
of the form 

e^v c (i-L/2) 2 , (7) 



with v c — 0.008. Choosing the parameter this way al- 
lows for nice fillings in the center and for densities going 
smoothly to zero near the edges of the trap. 

In Fig. [7| we see how plateaus with local fillings of an 
integer number of bosons can arise as more and more 
particles enter the system. The global compressibility is 
never zero, as it is the case in the unconfined model in 
the thermodynamic limit. We cannot speak therefore of 
a true quantum phase transition, the confining potential 
effectively cuts off all long range correlations. However, 
in local regions the local compressibility can get very 
low and the local density can get stuck at integer values, 
reflecting a local Mott region. This can be seen in 
Fig. 121 All these results are completely in line with 
those of Ref. |23|. Also, for a canonical calculation 
with an incommensurate filling a Mott phase with 
integer density can still be found, because the confining 
potential changes the local chemical potential. 

Looking along the sites can be interpreted as different 
/^-slices of the phase diagram in the (U, fi) plane for the 
unconfined model |(J . This allows to calculate the site at 
which a Mott domain is entered or left. It is also clear 
that the BKT transition has no analog in the unconfined 
case. We refer to Ref. [23j for a state diagram. The au- 
thors of Ref. [23 also claim that Ki ~ (rii — 1) as the Mott 
lobe is approached, independent of the on-site repulsion 
U or the chemical potential /i. In our calculations the 
same behavior was seen for parameters that are of the 
same order of magnitude, but for small and large values 
of U the local compressibilities did not reach to the same 
values in the Mott region. 
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FIG. 7: Profile of the local density n along the sites and as 
a function of the chemical potential fi. The on-site repulsion 
is U/2 — 7.1 Above a certain value of /x we see the emer- 
gence of a plateau (n = 1), and when the total number of 
particles is even more increased, we see the re-emergence of 
a compressible region. This happens first around the center 
and continues to exist till a plateau with n — 2 is reached. 
These results confirm the result of Ref. (2^1 
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FIG. 8: Profile of the local density n (solid line) and lo- 
cal compressibility k (dashed line) along the sites for U/2 = 
7.1, n/2 = 6.1, L = 128,« c = 0.008, confirming again the re- 
sult of Ref. 



B. Results in two dimensions 

1. Homogeneous case 

In principle it is possible to determine the phase 
diagram, but a complication that makes a comparison 
more difficult is that in the literature calculations 
are usually based on a fixed density while we are work- 
ing in the grand-canonical ensemble. The physically 
most interesting case is with a disordered chemical 
potential A phase diagram requires a study of 

phase transitions and very close to a transition point 



it is important to include more and more states into 
the truncated Hilbert space and at the same time going 
to larger lattices. The method described here can only 
give a qualitative answer and is not fit to quantitatively 
yield the exact location of the transition point and 
does not allow to calculate the critical exponents in an 
unambiguous way. 

The problem encountered here is a 'memory effect' 
when N s is not high enough. When the source terms are 
set to zero and a large lattice of L — 256 x 256 is taken, a 
calculation with a too low N s predicts a Mott phase while 
an increased N s leads to a SF phase. So, starting from a 
Mott phase (zero source terms), results in a Mott phase 
and starting form a SF phase (finite source terms) re- 
veals a SF phase. The issue of the phase diagram is very 
similar to the difficulties encountered with the strong- 
coupling expansion by Freericks et al. [fj , although their 
starting point is entirely different. As they point out, 
their method cannot describe the physics close to the 
tricritical point, the density fluctuations dominate even 
close to the tricritical point, and they can notice that the 
shape of the Mott lobes has changed from one to higher 
dimensions. Due to the limitations in our method we see 
the same qualitative aspects, but we ran into the same 
quantitative difficulties, with the same order of uncer- 
tainties. We will not report on calculations of the phase 
diagram here. 

2. Confined case 

The trapping potential takes on each site i the value 
£i = VcT^, (8) 

where measures the distance from the present site i to 
the center of the trap. The same holds in three dimen- 
sions. In Fig. [§] we plot the local density for a system 
of L = 64 x 64 sites, with U = 23.2,// = 28.0, v c = 0.05 
and the space is constantly truncated to 32 states, princi- 
pally in line with our philosophy of a limited but fast and 
reliable calculation. In Fig. ^| we show the local com- 
pressibility for the same system, but only one quarter of 
the figure is shown. The other parts are symmetric. Note 
again that there exists a Mott insulating region with in- 
teger density. The transition from the SF region to this 
Mott region is not sharp, and the local compressibility in 
the Mott region is small but remains finite. 

Another example can be found in Fig. 1111 and Fig. 1121 
for a lattice of 128 x 128 sites, showing Mott behavior 
and where for a slightly weaker U a new SF region would 
emerge in the center of the trap. 

C. Results in three dimensions 

The original experiments by Greiner et al. Q were 
performed in three dimensions, with laser beams cutting 
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FIG. 9: Local density for a lattice consisting of L = 64 x 64 
sites and with parameters U = 23.2, n — 28.0, v c = 0.05. Note 
again the region with fixed integer density and the smooth 
transitions. 




FIG. 10: Local compressibility as a function of the site indices 
and with the same parameters as in Fig. El 



the atomic cloud in about L = 65 x 65 x 65 sites, 
and a local density varying around n = 2.5 atoms per 
site. Mott and SF behavior were demonstrated after 
examination of the interference pattern of the laser 
images of the free expanding cloud. This means that 
the quantity of computational interest is the momentum 
distribution. In the original experimental setup, the 
absorption images of the three dimensional distribution 
are taken along two orthogonal axes, revealing only the 
integral over the third direction. The observed fading of 
the Bragg peaks had nothing to do with the appearance 
of Mott behavior, and happened actually when the 
system was already very deep in the Mott insulating 
phase. 

So what could be a clear signal of the transition? It 
was argued in Ref. JUJ that satellite peaks in the mo- 
mentum distribution are related to the appearance of a 




FIG. 11: Local density as a function of the site indices 
for a lattice of L = 128 x 128 sites and with parameters 
U — 22.0, /i = 35.6, v c = 0.008. The system is very close 
to developing a new SF peak in the center. 




FIG. 12: Same as in Fig. II II but for the local compressibility 
now. The parts of the plot that are not shown are symmetric. 



Mott region in the center of the trap. Once the Mott 
region spanned almost the entire lattice, the peaks disap- 
peared into the typical broad, low-peaked Mott distribu- 
tion. However, their Worm Monte-Carlo calculation was 
only on a lattice of L = 16x16x16 and it can be expected 
that for a larger lattice the central SF peak might be so 
dominant that the satellite peaks can hardly be resolved. 
We present a calculation on a lattice of L = 32 x 32 x 32. 
We show the momentum distribution in Fig. 1131 along 
the (1, 0, 0) axis, for a system with an emerging Mott re- 
gion of n — 1 in the center of the trap. As in Ref. [24j 
we see satellite peaks along the (1,0,0) direction, but 
the central peak dominates. The satellite peaks are only 
4.5% in magnitude of the central peak and will be dif- 
ficult to resolve in practice. For the experimental setup 
with its lattice of about L = 65 x 65 x 65 sites, the situ- 
ation will even be worse. We also note that the satellite 
peaks depend on the direction of investigation, no satel- 
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FIG. 13: Momentum distribution for a system with a Mott 
plateau (n = 1) in the center. The parameters are U/z = 
6.5, n/z = 2.6, v c = 0.04 and L = 32 x 32 x 32. The dashed line 
represents the distribution along the (1,0,0) direction, while 
the full line is taken along the (1,1,1) direction. According to 
Ref. [24| the satellite peaks in the dashed curve point at an 
emerging Mott region in the center of the trap. 



FIG. 14: Momentum distribution along the (0, 0, 1) axis for 
a system with particle densities varying between n = 2 and 
n = 3.2 (U/z = 11, At/2 = 30, v c = 0.1, L = 32 x 32 x 32). 
The appearance of satellite peaks cannot be related to the 
emergence of a Mott region in the center of the trap. 



0.68 




lite peaks were seen e.g. along the (1, 1, 1) direction in 
Fig. ^5] This direction dependency is a consequence of 
the breaking of rotational symmetry in a finite lattice, 
and its effects should diminish when larger lattices are 
taken. 

Furthermore, the average density in the experiments 
was about n w 2.5 in the center of the trap There 
are no satellite peaks when the central density is non- 
integer, despite a broad Mott n = 1 region for an on-site 
repulsion U that is strong enough. This Mott region is 
reflected in the tail of the momentum distribution |2 1| . 
For a system with local densities varying between n = 2 
and n — 3.2 (all non-integer densities), we nevertheless 
found satellite peaks in Fig. O and calculations showed 
the same behavior for densities ranging between n = 3 
and n = 4. These peaks cannot possibly be related to 
the emergence of a Mott region in the center of the trap. 
Local densities of n = 2 at the border of the trap cannot 
occur experimentally, but this situation can be thought 
of as the central region of a larger lattice, from which the 
outer regions are not trapped any more. 

When going to higher values of U and it is in 
principle possible to have Mott phases at n = 2 and 
n = 3. In the mean- field phase diagram of the homoge- 
neous model the different Mott lobes corresponding 
to densities n = l,n = 2, etc. get closer to each 
other along the direction of the chemical potential fi 
(see Fig. P). With the confining potential a present, 
the local densities along the different sites can be 
interpreted as a scan of the homogeneous model @- 
Hence in a small finite lattice it is not a priori clear 
if there are non-integer densities between the different 
broad Mott regions. For a system with parameters 
U/z = 30, fi/z = 75, v c = 9.1 and L = 16 x 16 x 16, we 
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FIG. 15: Momentum distribution along the (0, 0, 1) axis for 
a system with particle densities varying between n — at 
the edge of the trap to n — 3 at the center (U/z = 30, fj./z = 
75, v c = 9.1, L = 16 x 16 x 16), leading to 2280 particles in the 
system. The distribution is broad and not peaked, signaling 
virtually overall Mott behavior. 



found very few non-integer densities. The density profile 
consisted of four plateaus with n = 0, 1, 2, 3 respectively. 
We have almost a superposition of four Mott phases 
leading to the momentum distribution in Fig. 1151 which 
is very low peaked and broad. When the local particle 
density in the center of the trap is gradually increased 
from n < 3 till the Mott region with n = 3, and while 
there already exist broad Mott regions with n = 1 and 
n = 2, we did not witness any satellite peaks, because 
the Mott behavior of the n = 1 and n = 2 plateaus 
already dominated the momentum distribution. 

We are led to the observation that it will be difficult 
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to indicate the transition experimentally by satellite 
peaks, and that only examination of the intensity and 
the width of the central peak along one direction might 
be at hand to directly reveal Mott behavior. 



V. CONCLUSION 

In summary, we studied the Bose-Hubbard model 
subject to a confining field in the grand-canonical 
ensemble. We combined the Gutzwiller mean-field (MF) 
ansatz with a numerical (block) renormalization group 
method (NRG) and we could calculate observables like 
local densities, energies, the momentum distribution 
etc. The goal was to achieve variational results with 
energies much lower than in mean-field theory and at 
a low computational cost in order to make studies of 
large lattices in higher dimensions feasible. We have 
extensively discussed the advantages and limitations of 
this new method. The inclusion of source terms on the 
edges of the blocks improved results in the SF phase. 



We have examined the smooth transition between SF 
and Mott regions in the presence of a confining field. 
Although there is no real 'order parameter' to be found in 
the momentum distribution, the momentum distribution 
can nevertheless reveal important qualitative differences 
between pure SF systems and systems with dominant 
Mott behavior. These differences can experimentally best 
be seen in the central peak. Possible satellite peaks might 
be difficult to resolve when the total number of confined 
particles is large and when the filling factors are not of 
order unity. 
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